#ifndef RSIM_ODE23_H
#define RSIM_ODE23_H

/* -------------------------------------------------------------------------- *
 * File: Ode23.h                                                              *
 * Authors: Kishor Bhalerao                                                   *
 * Email : kishor8dm@gmail.com                                                *
 * Contributors:                                                              *
 * Email:                                                                     *
 * copyright (c) 2010 Authors.                                                *
 *                                                                            *
 * This program is free software: you can redistribute it and/or modify it    *
 * under the terms of the GNU General Public License as published by the Free *
 * Software Foundation, either version 3 of the License, or any later version.*
 *                                                                            *
 * This program is distributed in the hope that it will be useful, but        *
 * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY *
 * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License   *
 * for more details. (http://www.gnu.org/licenses/)                           *
 * -------------------------------------------------------------------------- */

#include "Integrator.h"
#include "Vector.h"
#include "LinAlgDecl.h"

#include <iostream>

using namespace std;

namespace RSIM{

/** This class largely mimics Matlab ode23. YDot is evaluated 3 times per time-step */
class Ode23:public Integrator{
        public:
                Ode23():Integrator(){}
                
                void initIntegrator(const Vector& Y);
                
                /** Error order of this integrator*/
                const int getErrorOrder()const{return 3;}
                        
                void attemptSingleStep(void(*calcYDot)(const double& ,	// Current time (t)
                                                const Vector&, 	// Vector Y(t)
                                                Vector& 		// Vector dY/dt
                                                ),
                                        const double& CurrentTime, // given: t
                                        const Vector& Y0, 	  // given: Y(t)
                                        const Vector& YDot,	  // given: YDot(t)
                                        const double& deltaT, 	  // given: deltaT
                                        Vector& Y_deltaT, 	  // calculate:Y(t+deltaT)
                                        Vector& YDot_TPlusDeltaT,  // calculate YDot(t+deltat)
                                        double& error		  // calculate: estimate of local error in this step.
                );
                
                void interpolate(const double& t0, 
                                const Vector& y0,
                                const Vector& y0Dot,
                                const Vector& yDot,
                                const double& h, 
                                const double& t_interp, 
                                Vector& y_interp);
                
        private:
                /** The space for k is allocated by stepTo function calling initIntegrator. */
                Vector k[2];
};// class Ode23

}// namespace RSIM

#endif